#ifndef RG_COREMATH_H
#define RG_COREMATH_H

#include "config.h"
#include <cmath>
using namespace std;

namespace rg{
  const real pi=acos(-1.0f);
  const real inf=1.0/0.0;
  const real eps=1e-12;

  struct vec3 {
    union{
      struct{real x,y,z;};
      struct{real elem[3];};
    };

    explicit vec3(const real& a=0,const real& b=0,const real& c=0):
      x(a),y(b),z(c){}

    real operator[](int i)const{return elem[i];}
    real& operator[](int i){return elem[i];}

    vec3 operator+=(const vec3& v){x+=v.x; y+=v.y; z+=v.z;  return *this;}
    vec3 operator-=(const vec3& v){x-=v.x; y-=v.y; z-=v.z;  return *this;}
    vec3 operator*=(const real& s){x*=s; y*=s; z*=s;  return *this;}
    vec3 operator/=(const real& s){x/=s; y/=s; z/=s;  return *this;}

    real sqlen() const {return x*x + y*y + z*z;}
    real len() const{return sqrt(sqlen());}
    vec3 norm() const{vec3 u=*this; u/=len(); return u;}
    vec3 normalize() {real l=len(); x/=l; y/=l; z/=l; return *this;}

    vec3 operator-() const{return vec3(-x,-y,-z);}
    vec3 operator+(const vec3& v)const {return vec3(x+v.x,y+v.y,z+v.z);}
    vec3 operator-(const vec3& v)const {return vec3(x-v.x,y-v.y,z-v.z);}
    vec3 operator*(const real& s)const {return vec3(x*s,y*s,z*s);}
    vec3 operator*(const vec3& v)const {return vec3(x*v.x,y*v.y,z*v.z);}
    vec3 operator/(const real& s)const {return vec3(x/s,y/s,z/s);}
  };
  inline vec3 operator*(const real& s,const vec3& v){
    return vec3(s*v.x,s*v.y,s*v.z);
  }
  inline vec3 cross(const vec3& u,const vec3& v){
    return vec3(u.y * v.z - u.z * v.y,
		u.z * v.x - u.x * v.z,
		u.x * v.y - u.y * v.x);
  }
  inline real dot(const vec3& u,const vec3& v){
    return u.x * v.x + u.y * v.y + u.z * v.z;
  }


  struct vec4{
    union{
      struct{real x,y,z,w;};
      struct{real elem[4];};
    };

    explicit vec4(const real& a,const real& b,const real& c,const real& d=1):
      x(a),y(b),z(c),w(d){}
    explicit vec4(const vec3& u=vec3(0,0,0),const real& s=1):
      x(u.x),y(u.y),z(u.z),w(s){}

    vec3 v()const{return vec3(x,y,z);}
    real operator[](int i)const{return elem[i];}
    real& operator[](int i){return elem[i];}

    vec4 operator-()const{return vec4(-x,-y,-z,-w);}
    vec4 operator+=(const vec4& u){x+=u.x; y+=u.y; z+=u.z; w+=u.w; return *this;}
    vec4 operator-=(const vec4& u){x-=u.x; y-=u.y; z-=u.z; w-=u.w; return *this;}
    vec4 operator*=(const real& s){x*=s; y*=s; z*=s; w*=s; return *this;}
    vec4 operator/=(const real& s){x/=s; y/=s; z/=s; w/=s; return *this;}

    real sqlen() const {return x*x + y*y + z*z + w*w;}
    real len() const{return sqrt(sqlen());}
    vec4 norm() const{vec4 u=*this; u/=len(); return u;}
    vec4 normalize() {real l=len(); x/=l; y/=l; z/=l; w/=l; return *this;}

    vec4 operator+(const vec4& u)const {return vec4(x+u.x,y+u.y,z+u.z,w+u.w);}
    vec4 operator-(const vec4& u)const {return vec4(x-u.x,y-u.y,z-u.z,w-u.w);}
    vec4 operator*(const real& s)const {return vec4(x*s,y*s,z*s,w*s);}
    vec4 operator*(const vec4& v)const {return vec4(x*v.x,y*v.y,z*v.z,w*v.w);}
    vec4 operator/(const real& s)const {return vec4(x/s,y/s,z/s,w/s);}
  };
  inline vec4 operator*(const real& s,const vec4& u){
    return vec4(s*u.x,s*u.y,s*u.z,s*u.w);
  }

}//namespace rg

#endif
